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ABSTRACT 

The convective period leading up to a Type la supernova (SN la) explosion is 
characterized by very low Mach number flows, requiring hydrodynamical meth- 
ods well-suited to long-time integration. We continue the development of the low 
Mach number equation set for stellar scale flows by incorporating the effects of 
heat release due to external sources. Low Mach number hydrodynamics equations 
with a time-dependent background state are derived, and a numerical method 
based on the approximate projection formalism is presented. We demonstrate 
through validation with a fully compressible hydrodynamics code that this low 
Mach number model accurately captures the expansion of the stellar atmosphere 
as well as the local dynamics due to external heat sources. This algorithm pro- 
vides the basis for an efficient simulation tool for studying the ignition of SNe la. 

Subject headings: supernovae: general — white dwarfs — hydrodynamics — 
nuclear reactions, nucleosynthesis, abundances — convection — methods: nu- 
merical 

1. Introduction 

Modeling the period of convection leading up to the ignition of Type la supernovae 
(SNe la) is critical to determining the distribution of hot spots that seed the subsequent 
explosion. Multidimensional simulations of SNe la explosions presently seed one or more 
hot spots at or near the center of a white dwarf and use a flame model to describe the 
subsequent evolution (see for example Ropke & Hillebrandt 2005; Gamezo et al. 2005; Plewa 
et al. 2004). Variations in the size, number, and distribution of these seeds can lead to large 
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differences in the explosion outcome (Niemeyer et al. 1996; Garci'a-Senz & Bravo 2005; Livne 
et al. 2005). The convective flows leading up to ignition have Mach numbers of 0.01 or less, 
with temperature perturbations of only a few percent (Woosley 2001; Woosley et al. 2004) — 
conditions that are extremely challenging for fully compressible codes. One-dimensional 
statistical methods (Wunsch & Woosley 2004) predict that off-center ignition is likely, but 
they cannot give information about the distribution of the hot spots. Only recently has 
progress been made in multidimensional modeling of convection in white dwarfs. The flrst 
such calculations (Hoflich & Stein 2002) evolved a two-dimensional wedge of the star for a few 
hours using an implicit hydrodynamics algorithm. Convective velocities of approximately 
100 km s~^ developed. They observed compression near the center of the star leading to 
slightly off-center ignition. Three-dimensional anelastic calculations (Kuhlen et al. 2006), 
showed a large-scale dipole flow dominating the evolution, leading to an off-center ignition. 
These simulations used a spectral decomposition, with a small portion of the center of the 
star removed due to the coordinate singularity at r = 0. Neither of these calculations 
operated at a Reynolds number large enough to see fully developed turbulence. Further 
three-dimensional studies are needed to see how robust this dipole flow is to rotation and 
convection at higher Reynolds and Rayleigh numbers. 

The goal of the present work is the development of a new multidimensional hydrodynam- 
ics algorithm capable of evolving the full star from the convective phase, though ignition and 
into the early stages of flame propagation. Long time integration is critical. As we showed 
previously (Almgren et al. 2006 — henceforth paper I), the low Mach number hydrodynamics 
equations provide an accurate description of flows with Mach numbers less than 0.2. By 
filtering out sound waves, the low Mach number approximation allows for much larger time 
steps (~ 1/M larger) than corresponding compressible codes. In contrast to the anelastic 
equation set, the low Mach number equations are capable of modeling flows with flnite- 
amplitude density and temperature perturbations. The only restriction is that the pressure 
perturbation be small. Furthermore, because the compressibility effects due to both the 
background stratiflcation and local heat release are included, the low Mach number equa- 
tions set can self-consistently evolve the expansion of a hydrostatic atmosphere due to heat 
release (Almgren 2000). 

In this paper, we continue the development of the low Mach number hydrodynamics 
algorithm to include the effects of heat release. An energy equation is added, and an earlier 
assumption from paper I is relaxed, now allowing the background state to vary in time. 
In § 2 we develop the low Mach number equation set. In § 3 the numerical methodology is 
explained. Comparisons to fully compressible calculations are provided in § 4 to demonstrate 
the accuracy and utility of our new algorithm. We conclude in § 5. 
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2. Low Mach Number Hydrodynamics 

In paper I, we derived a system of low Macli number equations for stellar atmospheres 
for which there was a time-independent background state. The necessary assumption for 
validity of this system was that the Mach number (M) of the flow be small. In this case, 
any pressure deviations from the base state pressure, which are O(M^), are also small. The 
perturbations of density and temperature need not be small. 

This system is valid for many low Mach number terrestrial and stellar flows, but fails to 
capture the correct atmospheric response to large-scale heating that radially shifts the entire 
atmosphere at and above the level at which the heating occurs. This was shown analytically 
for the terrestrial atmosphere, using the pseudo-incompressible approximation, in Bannon 
(1996). In Almgren (2000) it was shown that when time variation of the background state is 
correctly included, the solution calculated by the low Mach number equation set is identical 
to that reached by the fully compressible equation set. 

Physically, this is consistent with the interpretation of the low Mach number equations 
as representing instantaneous acoustic equilibration. If heating of a local parcel of fluid 
results in a large temperature perturbation from the ambient, the effect of the resultant 
acoustic waves is to return the parcel to pressure equilibrium with the fluid around it by 
expansion of the parcel. The density and temperature variations of the parcel relative to the 
ambient values may be large, but the parcel will remain close to pressure equilibrium. 

By contrast, when an entire layer of the atmosphere is heated, the acoustic equilibration 
process brings the entire layer to a new hydrostatic equilibrium. Consider a horizontally 
uniform atmosphere with heating uniformly applied throughout a horizontal layer. The 
response of the atmosphere will itself be horizontally uniform, and for positive heating each 
parcel of fluid above the heated layer will rise in its respective radial column. In equilibrium, 
given the assumption that no fluid is lost at the top boundary, then following an upward shift 
of each parcel, the mass of fluid above any given parcel will not have changed. If gravitational 
acceleration is effectively constant over the length scale of the base state displacement, 
then the weight of fluid above a parcel will not have changed. Thus for an atmosphere 
in hydrostatic equilibrium the pressure of each parcel will not have changed, although that 
parcel will have changed its radial location. In other words, the material derivative, rather 
than the time derivative, of the base state pressure must be zero. Numerical examples in § 4 
of this paper will confirm the time-dependent response of the base state as well as of the full 
state is necessary to correctly capture the atmospheric response to large-scale heat release. 

We note that the assumption of constant gravity over the displacement distance of the 
base state does not mean one must assume constant gravity over the full domain. In practice 
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the length scale of the base state adjustment is much smaller than the length scale of more 
localized motions. 

We now generalize the low Mach number equation set from paper I to allow for time 
dependence of the base state. We recall from paper I the fully compressible equations of 
motion in a stellar environment, but here we add an external heat source, -ffext, and again 
neglect compositional and reaction terms: 

| + V.(pU) ^ 0, 

- V ■ (pUU) + Vp = -pge, , (1) 



dt 
d{pE) 



dt 

and an equation of state 



V ■ {pVE + pV) = V ■ (kVT) - pg{V ■ e,) + pH,^, , (2) 



P = P{P,T) 



Here p, U, T, and p are the density, velocity, temperature and pressure, respectively, E = 
e + U-U/2 is the total energy with e representing the internal energy, and g{r) is the radially 
dependent gravitational acceleration (resulting from spherically symmetric self-gravity), 
is the unit vector in the radial direction, and k, is the thermal conductivity. The Reynolds 
number of flows in a typical white dwarf is sufficiently large that we neglect viscosity here, 
though viscous terms could easily be included in the model and the numerical methodology. 

Again we choose to work with enthalpy, h = e + p/p, rather than energy, replacing 
Eq. [2] above by 

^?|i,V.,U..)^^.p.. (3) 

where pH = pH^^x, + V ■ (k VT) represents the enthalpy source terms, and D / Dt = 9t + U ■ V 
represents the Lagrangian (or material) derivative. (For the purposes of this paper we could 
alternatively use the entropy equation, 

P^Dt = ' 

but for future work in which the source terms due to reactions are essential we will prefer 
the enthalpy formulation.) 

We recall from the low Mach number asymptotics in paper I that the assumption that 
M = |U|/c <C 1 is sufficient to decompose the pressure into a base state pressure, po, and a 
perturbational pressure, vr, i.e. 



p(x,r,t) = po(^i) + 7i"(x,r,t) , 
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where vr/po = O(M^). Here x represents the horizontal coordinate directions and r represents 
the radial direction. We define the base state density, po? by assuming hydrostatic equilibrium 
of the base state; this allows us to rewrite Eq. [1] in the form 

d{p\J) 



dt 

with no loss of generality. 



V-(pUU) + V7r = -(p-po)'7e, 



Continuing to follow the derivations of paper I but with a time-dependent base state, 
we rewrite conservation of mass as an expression for the divergence of velocity: 

V.U^-i^. ,4) 

Differentiating the equation of state, p = p{p,T), along particle paths, we can write 

Dp 1 fDp DT\ 



Dt pp \ Dt Dt 

with Pp = dp/dp\rp, and pt = dp/dT\^. 

An expression for DT / Dt can be found by applying the chain rule to the enthalpy 
equation (Eq. [3]): 

' ^'1-pK)^-^ph) , (6) 



Dt pcp V Dt 
where Cp = dh/dT\^ is the specific heat at constant pressure, and hp = dh/dp\j. for conve- 
nience. Substituting Eq. [6] into Eq. [5] and the resulting expression into Eq. [4] yields 

V.u = ^(^(i-pM-i):^ + ^(«!I^) . 

ppp \pcp J Dt ppp \ pCp J 

still with no loss of generality. Now, replacing p by po{r, t) and recalling the definition of H, 
we write the divergence constraint as 

V-U + «f% + U-Vpo) = — f^(V-(KVT) + pffext)) , 
\dt J ppp \pcp J 

where 



p'^CpPp J Tipo 

and Fi = d{\ogp) / d{\og p)\^. We recall from paper I that V ■ U + aU ■ Vpo can be rewritten 
as l/(3oV ■ (/3oU) where 

A(n<)=/3(0.<)oxp(|j^|^*') . (7) 
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Thus we can write the constraint as 

V ■ (/5oU) = Po{S - a^) . (8) 

For the purposes of comparing the fundamental hydrodynamic behavior of this low Mach 
number model to the established compressible formulation, we will from now on neglect 
thermal conduction. Summarizing the low Mach number equation set for this specialized 
case, with the momentum equation re-written as an evolution equation for velocity, we have 



dp 

at 



-V ■ (pU) 



^ ^ _V.(pU/.) + ^ + pi/.., (9) 

f = -U.VU-lv.-(^.e. , 
at p p 

V.(m = ;3.(./^„-^|^) , (10) 

where we define, for convenience, a = PT/{pCpPp)- 

This system differs from that in paper I in that now po and po are unknowns as well as 
p, ph, U, and vr. The equation of state was used to derive the constraint thus to include it 
here would be redundant. When reactions and compositional effects are included in future 
work, evolution equations for species will be added to this system and reaction terms will be 
added to the enthalpy equation and divergence constraint, but for the hydrodynamical tests 
we present here this system is sufficient. 

We follow the approach used in Almgren (2000) to compute the time evolution of the 
base state, recalling from the beginning of this section that the pressure of each parcel 
remains unchanged during base state adjustment, i.e., 

^ = 0. (11) 
Dt ^ ' 

We first calculate the radial velocity field, denoted Wq, that adjusts the base state. We 
decompose the full velocity field, U, into w^er and the remaining velocity field, U, that 
governs the more local dynamics, i.e., 

U(x, r, t) = «;o(r, t) + U(x, r, t) , (12) 

and write Eq. [10] in terms of wq and U, 

V ■ (/?oti^oe.) + V ■ (/5oU) = /5o (^ai/ext - 1^^^) ■ (13) 
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We then integrate Eq. [13] over a horizontal slab Qh x {r — h,r + h) to obtain 

iPo ^^^^ 

Assuming solid wall or periodic boundary conditions on the horizontal boundaries, or that 
the horizontal velocity decays sufficiently as we reach the horizontal boundaries, we can 
simplify the volume integrals into area integrals over Qh, 

[ ({f3oWo) + {PolJ)-er)d^^ = [ pJaH,^,--^^)drd^, (15) 

To define the model we want to partition the velocity field so that the horizontal average of 
the radial fiux due to heating be entirely incorporated into Wq rather than U, namely 

U ■ (ix = , 

Using this relationship, and taking the limit as — Eq. [15] can be simplified to 

^f*""'' - ^„f(;^,-J_%) . (16) 



dr \ Fipo dt 



where we define aH = 1/Area(r2//) /^^^((TiJcxt) dx. 

We can further simplify Eq. [16] by expanding d{PoWo) / dr = Podwo/dr + wodPo/dr, 
exploiting Dpo/Dt = dpo/dt + Wodpo/dr = to replace dpo/dt, and recalling from the 
definition of /3o that (l//3o) OPo/dr = (l/Fipo) dpo/dr (see Appendix B of paper I for the 
derivation of Po-) Then 

dwp ^ — 
dr 

This system can be integrated by noting that if there exists a lower boundary at r = tq with 
zero normal velocity (such as the center of a star), then at any time, t, 

wo{r,t) = / '^{r',t) dr' . (17) 

J ro 

The base state pressure and density update follow from Eq. [11] and conservation of mass, 
respectively: 

(18) 
(19) 



dpo 


dpo 


at 


or 






dt 


dr 
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There are now two choices for defining the new base state enthalpy; these options are an- 
alytically equivalent but may differ numerically. The first is to use the equation of state: 
{ph)o = po^(PO)Po)- The second is to use Eq. [9]. In this second approach, we can exploit 
Dpo/Dt = 0, but must also correctly partition the heating term in the right hand side of 
Eq. [9]. The partitioning is given by the requirement that the base state continue to satisfy 
the equation of state. Then, given Dpo/Dt = 0, 

Dho dho 



Dt dT 
1 



DTq Pp^DpQ poCpPp—^ 

— — — = Cr, ( I — — — = (J ri 

Dt Pt Dt Pt 



-aH , 

recalling a = PT/{pCpPp) and letting ctq = <^{PoyPo)- Returning to conservation form, we can 
write 

d{ph)o d{wo{ph)o) ^ Po —rj 

— ^7 — = ^ \ crii . [20) 

ot Or ao 

Using the velocity decomposition (Eq. [12]) we can re- write the evolution equations for 
p and ph as 

I . -V.(pU)-^ (21) 
^ = -V.(p/,,U)-^ + *^ + ,//„, (22) 

where w = U-e^. We can also write these in perturbational form (with no loss of generality): 
do' 

-f- = _v.(p'(U + ti;oe,))-V-(poU) 



V ■ ((p/i)'(U + WoBr)) - V ■ ((p/i)oU) + + pH,,, - ^aH 

ot Or \ ctq 

using Eq. [19] and Eq. [20], where p' = p — po and (ph)' = (ph) — {ph)o. 
The evolution of the velocity field becomes 

<9U ~ ~ d\J ^dwo 1 (p-po) /oQ^ 

— - = -U ■ VU - Wo— w^Br - -Vtt ger , (23) 

ot or or p p 

and subtracting Eq. [16] from Eq. [10], the constraint equation for U becomes 

V ■ (/3oU) = Poi^rH)' . (24) 



where we define {crH)' = cri^cxt — {cfH). 

In summary, then, the evolution of the base state is described by Eq. [18] and Eq. [19], 
with Wo given by Eq. [17], and the evolution of the full state is given by Eq. [21], Eq. [22] 
and Eq. [23] with the divergence constraint given by Eq. [24]. 
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3. Numerical Methodology 

Our strategy for evolving the low Mach number system with a time-varying base state 
is a fractional step approach. In each time step we first update density and enthalpy as if the 
base state were time-independent, giving us predicted values that can be used to construct 
time-centered values in the right hand side of Eq. [17]. We then compute the evolution of 
the base state and recompute the updates to density and enthalpy, incorporating the base 
state adjustment. Finally, we update and project the velocity field to define the new values 
of velocity and pressure. The upwind methodology used to update all the state variables 
provides a robust discretization of the convective terms that avoids any stability restriction 
other than the CFL constraint, i.e. the time step scales linearly with the grid spacing and 
inversely with the maximum magnitude of the velocity in any one coordinate direction in 
the domain. 

All base state quantities as well as all state quantities other than the perturbational 
pressure, tt, are defined at cell centers and integer time levels. The perturbational pressure 
is defined at nodes and at half-times; similarly, the advective velocity and fluxes used for 
advective updates are defined at edges and half-times. In this section we replace U by U for 
convenience of notation. 

Initialization Specification of the initial value problem includes initial values for po,Po 
and ho (or Tq) as well as U, p and h (or T) at time t = 0, and a description of the boundary 
conditions, but the perturbational pressure is not initially prescribed. We calculate Po at 
t = using Eq. [7]. Given this initial (3q, we project the initial velocity field to ensure that 
it satisfies the divergence constraint at t = 0. Then initial iterations of the following steps 
(typically two are sufficient) are performed to calculate an approximation to the perturba- 
tional pressure at t = At/2. At the end of each initial iteration all variables other than vr 
are reset to their initial values. 

The following steps are components of the single time step taken to advance the solution 
from to 

Step 1 In this step we construct U^°^, a time-centered, second-order accurate, staggered- 
grid approximation to U at t^^^^^, using an unsplit second-order Godunov procedure (Colella 
1990). To do so we first predict U^°^'* using the cell-centered data at and the lagged 
pressure gradient from the interval centered at t""'^^^. The provisional field, U^°^'*, represents 
a normal velocity on cell edges analogous to a MAC-type staggered grid discretization of the 
Navier-Stokes equations (Harlow & Welch 1965). However, U^^^'* fails to satisfy the time- 
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centered divergence constraint (Eq. [24]). We apply a discrete projection by solving the 
elliptic equation 

on 

^MAC^^^MAC^MAC) ^ ^MAC^^nuADV,*) _ ((^i/)')" 

for 0^^"^^, where _D^^c ^presents a centered approximation to a cell-based divergence from 
edge-based velocities, and G^^^ represents a centered approximation to edge-based gradients 
from cell-centered data. The solution, (p^^^, is then used to define 

-jjADV _ -jjADV,* ^^MAC jMAC 

pn 

In the above equations, we average /3q and Pq from centers to edges, i.e., Po'j+y^ ~ ^f^iPo] + 
/3o ■+!), Pr+y,,, = l/2(pi:, + and .^y^ = l/2(p^^. + pl^^,). 

Step 2 We update p and ph as if t^o = and the base state were constant, i.e., we discretize 

1^ = -V ■ (p'U) - V ■ (poU) , 

^ = -V-{{phyV)-V-{{ph)oV)+w^ + pH,^, 

using the second-order advection methodology as in paper I with pHcxt treated as an explicit 
source term. The discretization takes the form 

p"+^'* = p" -At[V-(p'U^°^)]"^'''^-AtV-(p^U^°^) , 
(p/,)«+i.* = (p/i)--At[V-((p/i)'U^°^)]"'''^'-AtV-((p/i)o"U^°^) 



where = U^°^ ■ e,. and p"+'/2 = l/2(p" + p'^+i-*) in the construction of {pH^xtT^^'\ 

The details of the upwind construction of [V ■ (U^^^p')]"^'''' and [V ■ (U^°^(p/i)')] are 
given in Appendix A, where we consider the construction of [V ■ (Vs)]""^^^ for any edge-based 
vector field V and cell-centered quantity s. In this step V = U^°^. The terms, V- (pqU^^^) 
and V ■ ((p/i)qU'^°^), are defined differently, in that we do not upwind po or (p/i)o in this 
step, rather they are simply averaged onto edges as /?o and p" were averaged in Step 1. 

Step 3 We integrate Eq. [17] to determine wq on edges. 
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using the equation of state given p"+^/2 g^j^^ compute a. We then update the base state 
quantities, 

(Po)r = (Po)--|^((PoH>|-(/'o^o);!^) , 

The construction of Po^^y^ and Po"+^^ is described in Appendix A. After construction of the 

new base state we compute using Eq. [7], then set = l/2{/3^ + We use the 

equation of state here to calculate (p/i)o'^^ in order to keep the base state thermodynamically 
consistent. 



Step 4 In this step we repeat the update of p and ph, but in the prediction of the edge 
states here V = XJ^°^ + woGr. We also center the w dpo/dr term in time: 

p"+^ = pr + (p"-pS)-At[V.(p'(U^°^ + t.oe.))]"^'/^-AtV.(pSU^°^) 

(p/^r+i = (p/.)r^ + ((p/.r-(p/z)^')-At[V-((p/.y(U^°^ + t.oe.O)]"^'^^ 
-AtV ■ ((p/i)^U^°^) 

A^ ADv [ f dpoY^^ , f 9PoY\ , f u Po 



We use the perturbational form of these equations in order to ensure that numerically, if, as 
in the anelastic case, Po = po, and p"^ = Pq and H is horizontally uniform, then p*^"*"^ = Pg"^^, 
i.e. no perturbation to the base state density is introduced in a case where analytically there 
should be none. 



Step 5 We then update the velocity field, U" to U"^^'* by discretizing Eq. [23], 

n+ 

Or 



p"+y2 pn+y2 ^ ' 

with p"+y2 = y2(^p"-|-p"+i) and G a discretization of the gradient operator. The construction 
of [((U^°^ + woer) ■ V)U]"^'''' is described in Appendix A with V = U^°^ + Woer and s 
set to each component of U" individually. Finally, we impose the constraint (Eq. [24]) 
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by solving 

for nodal values of (f) where is the standard bilinear finite element approximation to 
V ■ {Po/p)\/ with p and /3o evaluated at t"'"''/^. (See Almgren et al. (1996) for a detailed 
discussion of this approximate projection; see Almgren et al. (2000) for a discussion of this 
particular form of the projection operand.) We determine the new-time velocity field from 



Un+l _ jjn+1,* _ /q^ _ Q^n^y2\ 

£," + 72 V / 



At 

? 

and the new time-centered perturbational pressure from 



4. Numerical Results 

We consider three numerical tests in this section, each studying the response of the 
atmosphere to prescribed external heating. For each test case, the initial conditions in 
the computational domain are specified in two parts. The lower portion of the domain is 
initialized with a one-dimensional hydrostatic white dwarf model up until the outer boundary 
of the white dwarf. This initialization is identical for the compressible and low Mach number 
models. The model is created by specifying a base density of 2.6 x 10^ g cm""^ and base 
temperature of 7 x 10® K and integrating the equation of hydrostatic equilibrium outward 
while constraining the model to be isentropic. The composition is held constant at 0.3 ^^C 
and 0.7 ^^O, and the gravitational acceleration is fixed at —1.5 x 10^° cm s~^. We use the 
stellar equation of state developed by Timmes & Swesty (2000). This procedure provides a 
reasonable approximation of the state of the white dwarf just before runaway. None of the 
methods described here require constant gravity, but it was assumed for simplicity in the 
comparisons. 

The upper portion of the domain represents the region beyond the outer boundary of the 
white dwarf, and different approximations are used there for the compressible and low Mach 
number models. For the compressible calculations, the integration proceeds radially outward 
until the density reaches a threshold value of 10~^ g cm~^. Throughout the integration we 
set a low temperature cutoff of 10^ K, to keep the temperature in the outer layers of the 
model reasonable. Once the density drops below its cutoff, the integration is stopped and the 
material above it is held at constant density and temperature. This buffer region is necessary 
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to allow for expansion of the star; otherwise, as the star expands the loss of mass through 
the upper domain boundary would change the base pressure (Glasner et al. 2005), impacting 
the dynamics throughout the domain. Finally, for the multidimensional test cases, we add a 
convectively stable layer below the atmosphere to prevent any motions generated from the 
heating from interfering with the lower boundary. Figure 1 shows the initial temperature, 
density, entropy, and adiabatic indices (Fi and 7e = p/(pe) + 1) as a function of height for 
the compressible background. 

For the low Mach number model applied to the second and third test cases, the density 
cutoff is set to 2.5 x 10^ g cm~^, approximately the value at which the temperature cutoff is 
applied for the compressible background. Once the density reaches this cutoff the density, 
temperature and pressure are held constant, equivalent to gravity being set to zero radially 
outward of that position. Because the base state density in the buffer region is signifi- 
cantly higher for the low Mach number calculations than for the compressible calculations, 
this buffer region serves to damp motions that reach it without impacting the hydrostatic 
equilibrium of regions toward the center. Since the time step for the entire calculation is 
determined by the largest velocity in the domain, this damping is essential for low Mach 
number calculations in order to avoid excessively large velocities above the cutoff that would 
dictate an excessively small time step. An additional approximation in the outer region is 
that we set /5o = po for po < 5 x 10^ g cm~^ in order to suppress spurious wave formation at 
the outer boundary of the star. 

In the first test, a layer of the star is heated for 5.0 s with a heating profile, 

H = Hoexp{-{r-rof/W^) , 

with To = 4x 10'' cm, W = 10^ cm, and Hq = Ix 10^^ erg g~^ s~^. This energy generation rate 
is quite a bit higher than we would expect during the smoldering phase of the convection 
leading up to an SN la (Woosley et al. 2004), but necessary to see a response with the 
compressible code on a reasonable timescale. It also provides a more stringent test of the 
hydrostatic adjustment than a lower energy generation rate would give. Because the initial 
conditions and the heating are both one-dimensional, we use a reduced one- dimensional form 
of the equations to solve the systems. We contrast three different systems of equations. 

For this one-dimensional test, we compare the low Mach number results to those pro- 
duced by the fully compressible PPM method (Colella & Woodward 1984), as implemented 
in the FLASH code (Fryxell et al. 2000). Two versions of the low Mach number algorithm 
are used for this test. The first is the low Mach number equation set with time- varying base 
state (as described in § 2). The second is a formulation of the low Mach number equations in 
which the base state is time-independent (equivalent to the equation set present in paper I). 
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Figure 2 shows the density, temperature and pressure for the three solutions at t = 
5 s. All simulations were run on a uniform grid with 768 zones spanning 2.5 x lO^cm. 
The compressible and low Mach number expanding background solutions show excellent 
agreement. In the region of heating, the temperature has increased enormously, with a 
corresponding density decrease. The amplitude of the density decrease is much smaller, due 
to the degenerate nature of the equation of state. The inset in the density plot shows the 
density adjustment on a linear scale, showing the decrease at the heating height and an 
increase in the density above this. Both equation sets reach the same solution in response 
to the heating. The differences above 1.7 x 10^ cm are due to the different treatments of 
the upper boundary, and are not significant to the atmospheric dynamics. As a result of the 
energy deposition, the model expands by almost 10'' cm, or about 5%. In contrast, when 
the background is not allowed to adjust (dashed line), the low Mach number model fails 
to capture the expansion. While the temperature increases in the region of the heating, 
there is no expansion in the material above the heating layer. This is consistent with the 
point made in Bannon (1996) that the pseudo-incompressible equation set does not give the 
correct solution in the terrestrial atmosphere when the base state is not allowed to vary in 
time, and the demonstration in Almgren (2000) that the correct solution is found when the 
base state does absorb the horizontally averaged heating. 

Figure 3 shows the difference in density before and after heating. Here we see more 
clearly the adjustment of the density structure as a result of the localized heating. There is 
a slight rise in the density in the compressible solution below the heated layer. This small 
error is due to the difficulty in constructing a lower hydrostatic boundary that allows sound 
waves to leave the domain. For the low Mach number case, the state remains unchanged 
below the heating layer, as it should. The compressible and expanding background low Mach 
number solutions show a large increase in the density above the heating layer. This is not 
present in the low Mach number case where the background was not allowed to expand. 

In even this one-dimensional simulation, the choice of boundary conditions at the top 
and bottom boundaries of the computational domain is critical for the compressible cal- 
culations. The fully compressible code generates sound waves as it struggles to keep the 
hydrostatic solution steady on the grid. The boundary conditions must allow these distur- 
bances to leave the domain or they will corrupt the solution. We use a hydrostatic lower 
boundary, integrating the pressure and density in the ghost cells using a fourth-order recon- 
struction of the pressure (Zingale et al. 2002), where the temperature is kept constant in the 
ghost cells and the velocities are given a zero gradient. This provides pressure support to 
the material while allowing sound waves to leave the domain. Further robustness is obtained 
by computing the hydrostatic structure in the boundary from the initial base density and 
temperature, and keeping this structure fixed in time. At the upper boundary, we set the 
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density and temperature to our cutoff values and allow the ffuid to move out of the domain 
(with a zero gradient) but set inward velocities to zero. This has the effect of keeping the 
fluid velocities in this region small, and therefore, they do not dictate the time step. 

The low Mach number calculations, by contrast, are not sensitive to mass exiting the 
top boundary, as the hydrostatic equilibrium is incorporated into the base state, which is 
independent of the total mass. A simple outflow boundary condition is used at the top (with 
any inflowing velocities set to zero), and a reflecting boundary condition is used at the bottom 
boundary. Additionally, the compressible calculations need to resolve the scale height of the 
atmosphere very well to suppress any ambient velocities generated by slight imbalances of 
the pressure gradient and gravitational force (see Zingale et al. 2002 for a discussion of this). 
This is not the case with the low Mach number method, so we expect our low Mach number 
solutions to numerically converge at a lower resolution than the compressible solutions. 

For the second and third tests, we consider fully two-dimensional heating profiles, in 
which both parts of the new low Mach number algorithm are fully exercised. For these 
cases we compare the new low Mach number formulations with the compressible solution. 
In addition to the PPM algorithm, these multidimensional tests are also run with an unsplit 
compressible algorithm (Colella 1990), adapted to a general EOS, and incorporated into the 
FLASH framework. This is the same implementation of the unsplit method we described in 
paper I. The computational domain for these tests is 2.5 x lO^cm by 3.5 x lO^cm, spanned 
by a uniform grid with 640 x 896 zones. Periodic boundary conditions are used on the sides 
of the domain. 

In the second test, we specify three local regions of heating, designed to mimic "hot 
spots," but no heated layer as in the first case. In this scenario, while both the horizontal 
average and the local deviations from the horizontal average are non-zero, the deviations are 
much larger than the average, so the dominant effect is the local rather than horizontally 
averaged atmospheric response to the heating. For the first two seconds the heating profile 
has the form 



(Tj, of the perturbations are given in Table 1. After the first two seconds H is set to zero, 
and we continue to follow the evolution until t = 4 s. Figure 4 shows a time sequence 
from t = ltot = 4sof temperature contours in a subset of the domain spanning from 
5 X lO^cm to 1.4x lO^cm high. The temperature is an independent variable in the compressible 
calculations, but in the low Mach number model, for these examples, we evolve density and 
compute temperature from the equation of state using density and po- 




(25) 



where di = -\/(! 



X — XiY + {r — and the amplitudes, a,, locations, (xj,rj), and widths. 
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At early times, the three methods agree very well; at later times they diverge slightly. 
The vertical speed of the bubbles appears greatest with the low Mach number methodology 
followed by the unsplit compressible formulation; the PPM generates the slowest bubbles. 
We note that the difference in height between the FLASH PPM and unsplit methods is 
comparable to the difference between the unsplit and the low Mach number bubble heights. 
The precise reason for these differences is not yet completely understood; however, they serve 
to underscore the sensitivity of these flows and the difficulties in simulating them accurately 
with either the compressible or low Mach number approach. 

Figure 5 shows a resolution study of this second test case for each of the three meth- 
ods. The general trend one observes is that the location of the bubble rises with increased 
resolution. However, we also notice that the low Mach number model appears to converge 
to a solution at a lower resolution than either of the compressible models. This is likely due 
to the fact that hydrostatic equilibrium is guaranteed in the low Mach number method by 
the base state, while the compressible methods need considerable resolution just to keep the 
background medium quiescent. We also notice the difficulty that the FLASH PPM method 
has at the higher resolution, evidenced by the strong oscillations in temperature. This was 
also observed and discussed in paper I. 

In the third test, we add the heated layer of the first test case to the three "hot spots" 
of the second case, resulting in a case for which the horizontal average of the heating is larger 
than the perturbation from the average. The heating profile has the form 

H = Ho |exp(-(r - ro) V^^') + f (1 + tanh(2 - d.M)) | (26) 

where vq = 7.5 x 10^ cm and the amplitudes, a^, locations, (xi,rj), and widths, cTj, of the 
perturbations are as in the previous case and are given in Table 1. We apply the heating 
source for 2 s. As in the one-dimensional uniform heating case, we place the heating layer 
a bit above the lower boundary, so as to minimize contamination of the solution from lower 
boundary effects. Figure 6 shows the temperature contours for the unsplit compressible and 
low Mach number solver at t = 1.5 s, t = 1.75 s and t = 2 s. We again notice excellent 
agreement between the low Mach number and compressible results. We do not expect the 
exact shape of the rising bubbles to match precisely given the extremely unstable nature of 
the bubble's surface, but there is overall agreement between methods. Once again the low 
Mach number bubble rises slightly faster. We do not show the FLASH PPM results here, 
as the noise resulting from the dimensional splitting dominates the solution. We also note 
that at the final time, the low Mach number result shows a disturbance in the upper left 
corner of the domain. This disturbance is a result of a spurious wave generated at the outer 
boundary, which is not present in the compressible results because of the different treatment 
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of the outer boundary condition. The question of how best to represent the outer boundary 
of the white dwarf in a low Mach number calculation is still an open research question; at this 
point we note that while the high-velocity disturbance does restrict the time step, it does not 
appear to impact the solution in the primary region of interest. We also expect that more 
realistic heating profiles in three-dimensional geometries will not generate the disturbances 
in the outer boundary that we have observed in this final test case. 

Figure 7 shows the horizontal average of the difference between density at t = 2s and 
t = Os for the low Mach number and unsplit compressible results shown in Figure 6. On the 
left of this plot we see that the compressible solution has an increase in density over time 
below the level of the applied heating; this was also present in the first test case, although to 
a lesser degree. This error results from the difficulty in prescribing an accurate hydrostatic 
lower boundary. We note, however, that this is less than a 0.1% relative error in the density 
at this boundary. At approximately the center of the heated layer (7.5 x 10^ cm) each 
calculation shows a negative average density variation. The small difference here between 
the compressible and low Mach number results here may also be a product of the lower 
boundary condition. Overall, however, the two algorithms agree well in the average response 
of the atmosphere to the heating. 

Finally, since computational efficiency as well as accuracy is necessary for successful long- 
time integration, we comment on the relative efficiency of the low Mach number algorithm. 
For the second test case, for example, to evolve the state to 4.0 s, the unsplit compressible 
method took 14272 time steps while the low Mach number algorithm took only 233 time 
steps. Evolving to just 2.0 s on a single processor (2.8 GHz Intel Xeon) using the Intel 9.0 
compilers (with -03 -ipo optimization flags) took 71.76 hours with the FLASH PPM solver. 
FLASH was setup to use 32 x 32 zone blocks for greatest efficiency. By comparison, the low 
Mach number method took only 0.46 hours — a factor of over 150 speed-up. The unsplit 
method takes the same timestep as the PPM algorithm but is approximately a factor of two 
slower, due to the additional work per timestep. 

For the final test case, evolving the solution to 1.5 s took 6532 timesteps for the unsplit 
method and 749 time steps for the low Mach number algorithm. Already by this point 
in the calculation, the spurious disturbance at the outer boundary of the star in the low 
Mach number algorithm is impacting the time step, decreasing the relative advantage of 
the low Mach number algorithm. As noted above, we expect this not to be the case for 
more realistic ignition scenarios. However, this case points out that to achieve the gains in 
efficiency possible with a low Mach number model one must successfully address this issue. 
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5. Conclusions 

We have introduced a new algorithm for evolving low Mach number flows in the presence 
of local and large-scale heating. By contrast with the previous low Mach number model, 
this new model allows time variation of the base state in order to account for atmospheric 
expansion due to large-scale heat sources. The time evolution of the base state must be 
calculated at each time step in addition to the local dynamics. Numerical comparisons of low 
Mach number simulations with simulations using a fully compressible code demonstrate that 
the low Mach number algorithm with a time-dependent base state can accurately capture 
the hydrostatic adjustment of an atmosphere as well as local dynamics in response to large- 
and small-scale heat release. 

Our long-term goal is to develop the capability for full star simulation using the new 
low Mach number approach. The fundamental low Mach number approach has been vali- 
dated with a number of simplified test cases, but further development is necessary to begin 
to perform detailed physics investigations of ignition and other problems of interest. This 
development will include extension to three dimensions with adaptive griding, radial repre- 
sentation of gravity and the base state within the three-dimensional setting, non-constant 
gravity, and the calculation of internal heating due to reaction networks. 

The tests presented above are quite demanding, and provided a challenge to both the 
low Mach number and the compressible solvers. The energy generation and resulting tem- 
perature/density contrasts during the convective phase of an SN la are much smaller. Based 
on the agreement demonstrated on these difficult tests, we are confident that the low Mach 
number hydrodynamics method will be a useful and efficient tool in exploring the problem 
of SNe la ignition. In addition, this algorithm is applicable to a wide range of problems 
outside of our target application (SNe la), including Type I X-ray bursts, classical novae, 
and convection in stars. 
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A. Construction of Advective Updates 

Consider tlie construction of an advective update in the form [V ■ [sV)]^^^', given the 
cell-centered velocity field U" = {u,v), an edge-based velocity field, V = (V^,V^) and a 
cell-centered scalar, s. For simplicity we will present the construction in two dimensions, 
although extension to three dimensions is straightforward and is given in detail in (Almgren 
et al. 1998). 

We first extrapolate s from cell centers at to edges at ^"+^2 ^ging a second-order 
Taylor series expansion in space and time. The time derivative is replaced using the evolution 
equation for s. If, for example, St = —V ■ (sV) = —V ■ Vs — sV ■ V, then 

-5,; I 1/ ~ Si j + ——Sx + —St 



extrapolated from («, j), and 

Ax At 

= Si+i^j + (-^ - + - Si+ij(V^ + v;)i+ij 

extrapolated from (i + 1, j). In evaluating these terms the first derivatives normal to the face 
(in this case sj^™) are evaluated using a monotonicity-limited fourth-order slope approxima- 
tion (Colella 1985). The construction of the transverse derivative terms (wsy in this case) 
are given in detail in (Almgren et al. 1998). Analogous formulae are used to predict values 
for s^^.^y and s^^^y at the other cell edges. 

Upwinding is used to determine s at each edge as follows: 

if V^^i . > 

%-\-— J 

and similarly for defining Sj^+i using V^. Finally, we define the conservative update term, 
[V ■ (sV)]-^ = (V^.,^S.,., - Vl.,s^.^„) + (yi,^kH^ - V^_.5.,_.) . 
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The construction of pn^.~^y and pn^-~^v in Step 3 is similar but not identical to the above 
procedure. Here we make no reference to U". Rather 

,Ar At 



extrapolated from j, and 



Ar At, 



= Po"+i - + ^oj+i— )(poDi+i 
extrapolated from j + 1. Here Wqj = l/2(woj+y2 + ""^Oj-i/a)- Upwinding then determines 

f Pof+y, if wo,+y, > 

[ PoJ+y, if w^oj+y^ < 

The evolution equation for po differs from that for po so the construction of Po^^yl differs 
slightly from that for Po'^lj^- 

R n / Ar At lim\ n/ \ 

POj+y, = POj + - ^Oi^)(por )j - -^Poj i^Oj+y, - woj^y,) 
extrapolated from j, and 

T n ,Ar \ / lim\ n / \ 

POj+y^ = POj+i - + w^0i+i^)(P0r. )i+i - - ^Oi+yJ 

extrapolated from j + 1. The upwinding procedure is the same. 
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Table 1: Location of heating sources for non- uniform heating terms, Eq. 25 and Eq. 26. 
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3x10' 



Fig. 1. — The white dwarf atmosphere initial modeL Shown are the density (top left), 
temperature (top right), entropy (bottom left), and adiabatic indices (bottom right). To 
prevent convective motions from hitting the lower boundary in our multidimensional tests, 
the first 5 x 10^ cm of this model is constructed to have a convectively stable entropy profile. 
The one-dimensional tests do not use this portion of the model. 
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Fig. 2. — Hydrostatic adjustment problem with uniform heating: density (top), temperature 
(middle) and pressure (bottom). The initial conditions are shown in gray. The solid black line 
is the fully compressible solution, the dotted line is the low Mach number formulation that 
allows for the base state expansion, and the dashed line is the low Mach number formulation 
assuming a fixed base state. The compressible and expanding base state low Mach number 
solutions show excellent agreement. The low Mach number model with a fixed base state is 
unable to capture the correct solution. The inset in the density plot shows the structure in 
the vicinity of the local heating. 
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Fig. 3. — Hydrostatic adjustment problem with uniform heating: the difference in density 
between t = and t = 5 s. The sohd black line is the fully compressible solution, the 
dotted line is the low Mach number formulation that allows for the base state expansion, 
and the dashed line is the low Mach number formulation assuming a fixed base state. We 
see close agreement between the compressible and low Mach number formulation with the 
time-dependent base state. 
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Fig. 4. — Second test case: temperature contours for the low Mach number (green), unsplit 
(red), and PPM (blue) solvers, shown at 1, 2, 3, and 4 s. Here, a heating source term 
gradually adds energy at three points in the domain (see Eq. 25) during the first 2 s of 
evolution. This gives rise to the three buoyant plumes seen in the panels. Contours span 
10^ K to 8 X 10® K, spaced every 5 x 10'^ K. 
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Fig. 5. — Resolution study for the second test case. Each pane shows the middle bubble from 
Figure 4 only; PPM results are in the left pane, unsplit results are in the middle pane, and 
low Mach number results are in the right pane. Within each pane, the left half is identical 
to the left half of the bubble in Figure 4, with the same color scheme. On the right half of 
each pane is a reflection around the center line of the bubble of the comparable results but 
at both a lower (320 x 448) and higher (1280 x 1792) resolution. The low resolution is shown 
in purple and the high resolution is shown in black. 
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time = 1 .5 s 
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Fig. 6. — Third test case: temperature contours for the low Mach number (green) and unspht 
(red) solvers. Here, in addition to the localized heating from the previous test (see Figure 4), 
there is a uniform heating layer centered at a height of 7.5 x 10^ cm, as specified by Eq. 26. 
Contours span 10^ K to 8 x 10^ K, spaced every 5 x 10'' K. 
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Fig. 7. — Third test case: horizontal average of the difference between density at t = 2 
and density at t = 0. The dotted hne shows the low Mach number results; the solid line 
shows results using the unsplit compressible formulation. The discrepancy between the two 
solutions on the left is due to the difficulties in prescribing an accurate hydrostatic lower 
boundary for the fully compressible calculation. Overall, however, the two algorithms agree 
well in the average response of the atmosphere to the heating. 



